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Unique features of non-magnetic insulator phase are revealed, and the phase diagram of the t — t' 
Hubbard model containing the diagonal transfers t' on a square lattice is presented. Using the path- 
integral renormalization group method, we find an antiferromagnetic phase for small next-nearest 
neighbor transfer t' and a stripe (or collinear) phase for large t' in the Mott insulating region of 
the strong onsite interaction U . For intermediate t' /t ~ 0.7 at large U /t > 7, we find a longer- 
period antiferromagnetic-insulator phase with 2x4 structure. In the Mott insulating region, we 
also find a quantum spin liquid (in other words, a non-magnetic insulator) phase near the Mott 
transition to paramagnetic metals for the t — t' Hubbard model on the square lattice as well as on 
the anisotropic triangular lattice. Correlated electrons often crystallize to the Mott insulator usually 
with some magnetic orders, whereas the "quantum spin liquid" has been a long-sought issue. We 
report numerical evidences that a nonmagnetic insulating (NMI) phase gets stabilized near the Mott 
transition with remarkable properties: The 2D Mott insulators on geometrically frustrated lattices 
contain a phase with gapless spin excitations and degeneracy of the ground state in the whole 
Brillouin zone of the total momentum. The obtained vanishing spin renormalization factor suggests 
that spin excitations do not propagate coherently in contrast to the conventional phases, where 
there exist either magnons in symmetry broken phases or particle-hole excitations in paramagnetic 
metals. It imposes a constraint on the possible pictures of quantum spin liquids and supports an 
interpretation for the existence of an unconventional quantum liquid. The present concept is useful 
in analyzing a variety of experimental results in frustrated magnets including organic BEDT-TTF 
compounds and ''He atoms adsorbed on graphite. 

PACS numbers: 71.30.+h, 71.20.Rv, 71.10.Fd, 75.10.Jm, 71.10.Hf 



I. INTRODUCTION 



Among various insulating states, those caused by elec- 
tronic Coulomb correlation effects, called the Mott insu- 
lator, show many remarkable phenomena such as high-Tc 
superconductivity and colossal magnetoresistance near 
it 0,l3]- However, it has also been an issue of long debate 
whether the Mott insulator has its own identity distin- 
guished from insulators like the band insulator. This is 
because the Mott insulator in most cases shows symme- 
try breakings such as antiferromagnetic order or dimer- 
ization, where the resultant folding of the Brillouin zone 
makes the band full and such insulators difficult to distin- 
guish from the band insulators because of the adiabatic 
continuity. 

Except for one-dimensional systems, the possibility of 
the inherent Mott insulator without conventional orders 
has been a long-sought challenge. The Mott insulator 
on the triangular lattice represented by the Heisenberg 
spin systems was proposed as a candidate ;3] . Although, 
the triangular Heisenberg system itself has been argued 
to show an antiferromagnetic (AF) order |j|, intensive 
studies on geometrical frustration effects have been stim- 
ulated. In particular, the existence of a quantum spin liq- 



uid phase has been established in recent unbiased numer- 
ical studies performed on two-dimensional lattices with 
geometrical frustration effects^ ^ JJ. The spin liquid 
has been interpreted to be stabilized by charge fluctua- 
tions enhanced near the Mott transition. 

Recently extensive experimental studies on frustrated 
quantum magnets such as those on triangular, Kagome, 
spinel and pyrochlore lattices H lol. ITol flTlIll ll3L 111 
as well as on triangular structure of "^He on graphite [l5j 
have been performed. They tend to show suppressions 
of magnetic orderings with large residual entropy with a 
gapless liquid feature for quasi 2D systems or "spin glass- 
like" behavior in 3D even for disorder-free compounds. 
These gapless and degenerate behaviors wait for a con- 
sistent theoretical understanding. 

In this paper, by extending and reexamining the pre- 
vious studies [1, la 13 > ^^ further show more detailed 
numerical evidence for the existence of a new type of in- 
herent Mott insulator near the Mott transition; singlet 
ground state with unusual degeneracy, namely gapless 
and dispersionless spin excitations. Our present results 
offer a useful underlying concept for the understanding 
of the puzzling feature in the experiments. 

In this paper, we also show that several different anti- 
ferromagnetic phases appear in the region of large U/t. 



This includes normal antiferromagnetic order stabilized 
for small t' /t with the Bragg wavenumber Q — (tTjTt), 
collinear (stripe) order for large t' /t with Q — (0,7r) 
and longer-period antiferromagnetic order for interme- 
diate t' /t with Q = (7r,7r/2). 



II. FRUSTRATED HUBBARD MODELS 

In the present study, we investigate the Hubbard model 
on two-dimensional frustrated lattices. The Hamiltonian 
by the standard notation reads 

H = Hk + Hu 
Hk ^ - J2 ^ (4c,. + H.c) + ^ t' (cLq. + H.c.) 

{ij),c {k,l),(T 

on a A''-site square lattice with a nearest neighbor (t) and 
two choices of diagonal next-nearest neighbor (t') transfer 
integrals in the configuration illustrated in Fig. ^^A] and 
[B]. The energy unit is taken by t. Hereafter we call the 
Hubbard models with the lattice structures illustrated in 
Fig. ^A] and [B] , the models [A] and [B] , respectively. 
The t' transfer integrals bring about geometrical frustra- 
tion. The i, j represent lattice points and cj^ {cja) is a 
creation (annihilation) operator of an electron with spin 
a on the i-th site. 



been introduced into the PIRG method and symmetries 
of the model can be handled explicitly and precision of 
solutions can be significantly enhanced |20l |. 

Here we briefly introduce the PIRG method and quan- 
tum number projections. The ground state [ipg) can be, 
in general, obtained by applying the projector e^"^^ to 
an arbitrary state | (^initial) which is not orthogonal to the 
true ground state as 



1^9 



-,-tH\ 



^initial/ 



(2) 



To operate exp[— ri?], we decompose exp[— T-ff] into 
exp[-riJ] ~ [exp[-ATiJA'] H* exp[-Ar7?[/J]-^ for smaU 
Ar, where r = AfAr. When we use the Slater 
determinant as the basis functions, the operation of 
exp[— AT_ffif] to a Slater determinant simply transforms 
to another single Slater determinant. On the other hand, 
the operation of exp[— Ari?c/.] can be performed by 
the Stratonovich- Hubbard transformation, where a single 
Slater determinant is transformed to a linear combination 
of two Slater determinants. Therefore, the operation of 
exp[— r_ff] increases the number of Slater determinants. 
To keep manageable number of Slater determinants in 
actual computation, we restrict number of Slater deter- 
minants by selecting variationally better ones. This pro- 
cess follows an idea of the renormalization group in the 
wavefunction form. Its detailed algorithm and procedure 
are found in Ref. [lal- After the operation of exp[— rif], 
the projected wave function can be given by an optimal 
form composed of L Slater determinants as 
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FIG. 1: Lattice structure of geometrically frustrated lattices 
on a square lattice. The nearest- and next-nearest-neighbor 
transfers are denoted by t and t' , respectively. 



III. PATH-INTEGRAL RENORMALIZATION 
GROUP METHOD 

The model requires accurate and unbiased theoretical 
calculations because of large fluctuation effects expected 
from the low dimensionality of space and the geometri- 
cal frustration effects due to nonzero t' . Recently, the 
path-integral renormalization group (PIRG) method [Ig 
opened a way of numerically studying the models with 
the frustration effects more thoroughly without the neg- 
ative sign problem and without relying on the Monte 
Carlo sampling. The efhciency of the method was es- 
tablished through a number of applications j^ |a, [13, [l3 • 
Furthermore, a quantum number projection method has 



(^h-^, 



1^^"^) 



aW\ 



(3) 



aW\ 



where c^'s are amplitudes of |(/)q ). Operation of the 

ground-state projection can give optimal c^'s and |(/)q )'s 
for a given L. 

In most cases, eq. (3) can give only upper-bound of 
the exact energy eigenvalue. Therefore, to obtain an ex- 
act energy, we consider an extrapolation method based 
on a relation between energy difference bE and energy 
variance A.E [la, [Q ■ Here the energy difference is de- 
fined as &E = {&) — {H)g and the energy variance is 



defined as AE 



H) 



Here, 



H)g stands for the 



true ground-state energy. For h/;^ '), we evaluate the 
energy £'(^) and energy variance AE'-^\ respectively. 

If IV''-^-') is a good approximation of the true state, 
the energy difference SE^^' is proportional to the en- 
ergy variance AE'-^\ Therefore extrapolating iJ^^^ into 
/\£;(i) ^. by increasing L systematically, we can accu- 
rately estimate the ground-state energy. 

Next we consider a simple combination of the PIRG 
and quantum number projection. The PIRG gives ap- 
proximate wave function for a given L which is composed 

of L linear combinations of (t>a )■ Though spontaneous 



symmetry breaking should not occur in finite size sys- 
tems, symmetries are sometimes broken in PIRG calcula- 
tions because of the limited number of basis functions, if 
symmetry projections are incompletely performed. How- 
ever, extrapolations to the thermodynamic limit recover 
the true ground state, in which possible symmetry break- 
ings are correctly evaluated. 

For finite size system, to handle wave functions with 
definite and exact symmetries, we can apply quantum 
number projection to this wave functions as 



V'(^)) = ^c„p|0i^) 



(4) 



where P is a quantum-number projection operator. We 
can use the same amplitudes Cq's and the same bases 



aW\' 



s which the PIRG determines, while this ampli- 



tude Cq 's can be easily reevaluated by diagonalization by 
using quantum-number projected bases, that is, we deter- 
mine Cq's by solving the generalized eigenvalue problem 
as 



Hap^ — ^ap^^ 



(5) 



where iV^^ 



(0^1 P |0„), Hip = {M HP |0„). The lat- 
ter procedure gives a lower energy eigenvalue. By adding 
this procedure for the PIRG basis, we evaluate the pro- 
jected energies and energy variances, -Eproj and AE^^.^- 
for each L. We can estimate accurate energy by extrap- 
olating the projected energy into zero variance. Conse- 
quently we can exactly treat the symmetry and extract 
the state with specified quantum numbers by the PIRG. 
We call this procedure PIRG+QP. 



IV. PHASE DIAGRAM ON t' - U PLANE 

For i' = 0, the metal-insulator transition occurs at 
U = 0. The t' offers additional dimension for pa- 
rameter spaces of Hubbard models. Recently the exis- 
tence of a nonmagnetic insulator (NMI) near the metal- 
insulator transition boundaries was reported [5|,|6| on the 
two-dimensional frustrated Hubbard model on a lattice 
by using the PIRG method. In the present paper, we 
thoroughly investigate two-dimensional parameter spaces 
spanned by U/t and t' /t. 

To obtain the ground state for each point ( U/t, t' /t ), 
we carry out calculations for 6x6, 8x8 and 10x10 lattices 
by the PIRG-f-QP method. By extrapolating the ground 
state energies per lattice site, we determine its thermo- 
dynamic value. By a thorough search of the parameter 
space and improved accuracy of the computation, the 
phase diagram is clarified for the model [A] in a wider re- 
gion of the parameter space than the previous study [3 . 
It becomes quantitatively more accurate and contains a 
new feature, which has not been revealed in the previous 
results 5]. In Fig. [2 we show the phase diagram, where 



various types of phase appear. For small U/t, paramag- 
netic metallic phase appears and for larger U/t and small 
t' /t, antiferromagnetic insulating phase appears. As t' /t 
increases, non-magnetic insulating phase appears. This 
feature has already been reported in Ref . p[ . In addition 
to these phases, we find two additional phases: one is 
another type of antiferromagnetic insulating phase and 
the other is stripe ordered insulating phase. The nature 
of these new phases will be reported in the following sec- 
tions. 

We determine the phase boundary in the following 
way. The ground state energies per site in the ther- 
modynamic limit are obtained for each mesh point of 
this parameter space at {U/t, t' /t) = (ni, n2 x 0.1) where 
ni,n2 = 1,...10 after the size extrapolation. The phase 
boundary is determined by interpolating these energies 
per site. 
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FIG. 2: (color online) Phase diagram of the Hubbard model 
with the lattice structure illustrated in Fig. 1 [A] in the pa- 
rameter space of \J scaled by i, and the frustration parameter 
t! /t. AFI, AFI2, Stripe, PM, and NMI represent two types of 
antiferromagnetic insulating, stripe shape insulating, param- 
agnetic metallic and nonmagnetic insulating phases, respec- 
tively. 



V. AFI PHASE 

For t'/t = 0, the AFI phase of the configuration illus- 
trated in Fig.|ni(A) appears. By the PIRG-I-QP method, 
we can evaluate an energy gap between the ground state 
{S — 0) and the first excited state {S = 1) for AFI phase. 

The finite-size gap for an ^ x ^ system in the chiral 



perturbation theory [2j| in the form 



AE- 



3.900265c 

Anpi 



Oi^)] 



(6) 



is fitted with the calculated results in Fig. 2| The finite- 
size gap nicely follows the form © in the AFI phase. For 
example, at U = 4, t = 1 and t' = 0, the fitting in Fig. 01 
shows the spin wave velocity c ~ 0.74 and the spin stiff- 
ness p ~ c/8.15, which are equivalent to the estimate of 
the Heisenberg model at the exchange coupling J = 0.45 
in the spin wave theory. The fitted values of c and J 
well reproduce the previous estimates (ex. J ^ 0.4) ob- 
tained from the susceptibility and the staggered magne- 
tizations 23]. For U/t — 6.0, we obtain the spin wave 
velocity c ^ 0.97 and the spin stiffness p ^ c/5AA, while 
for U/t = 8.0, c ^ 1.03 and the spin stiffness p ~ c/4.61. 
The excitations in the AFI phase well satisfy the tower 
structure of the low-energy excitation spectra based on 
the nonlinear sigma model description. The AFI phase 
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To examine this configuration, we consider the equal- 
time spin structure factor defined, in the momentum 
space, by 



Siq) 



N 



3N ^ ^^' ■ ^■" 



Jq{Ri-Rj) 



(7) 



where Si is the spin of the i-th site and Ri is the vector 
representing the coordinate of the i-th site. In Fig.l^B), 
momentum dependence of S{q) for the ground state at 
U/t = 9.0 and t' /t = 0.7 is shown, where we see distinct 
peaks at (tt, 7r/2) and (tt, 37r/2). In Fig. we plot finite- 
size scaling of these peak amplitudes, which indicates the 
existence of the long-ranged order. A usual AFI phase 
has a configuration pattern in Fig. |3fA), which has 2x2 
super structure. On the other hand, this new AFI phase 
has 2x4 super structure. Therefore, for smaller lattices, 
it is difficult to identify it. For An x 4n {n is an integer) 
lattice, this structure is naturally realized and peak po- 
sitions of spin correlation sharply appear at (tt, 7r/2) and 
(tt, 37r/2), while for (4n -f- 2) x (4n 4- 2) (n is an integer) 
lattice, peak positions of spin correlation become wider 
from (tt, 7r/2 — 5) to (tt, 7r/2 -|- 5) and from (tt, 37r/2 — 5) 
to (tt, 37r/2 4- 5) [6 is around 7r/4). Thus there is an ir- 
regularity in spin correlation as a function of lattice size. 
In Fig. IHl we show summed amplitude of 5(g) over the 
peak for 6x6 and 10 x 10 lattices. 



FIG. 3: Configurations of AFI phase (A), new AFI phase (B) 
and Stripe phase (C). 
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FIG. 4; (color online) Size scalings of the energy gaps for total 
spin 5=1,2 and 3 in the AFI phase (A) ([/ = 4, t = 1, t' = 0) 
and (B) ([/ = 6, t = 1, i' = 0). The solid curves are fitting by 
the form (6) and the dashed curves illustrate curves obtained 
from the S = Q fitting multiplied with the factor 45(5-1- 1)/3. 

is extended in the region of nonzero and moderate am- 
plitude of t' /t with large U/t in the phase diagram. 



VI. AFI PHASE WITH LONGER PERIOD 

In Fig. 2, AFI phase with longer period (AFI2) appears. 
Its schematic configuration is depicted in Fig.jSjB). 




FIG. 5: Equal-time spin correlations in the momentum space 
on 10 X 10 at half fiUing for AFI phase (A) ((/ = 6.0, i'/i = 
0.5), new AFI phase (B) {U/t = 9.0,*'/^ = 0.7), stripe 
phase (C) {U/t = 9.Q,t' /t = 1.0) and NMI phase (D) 
((7/t = 6.0,i7f = 0.7). 

For each basis |(/>) of Eq. Q, we apply shift op- 
eration Ta^Aj on the lattice and evaluate an overlap 
{4>\T/^i^^j\4>). This overlap becomes quite large when 
Ai = 2 or Aj = 4, which reflects a basic configuration in 
Fig. El B. 

In the strong-correlation limit {U/t — > oo), the frus- 
trated Hubbard models become the J1-J2 Heisenberg 
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Da,p = {OaOp) 
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FIG. 6: Finite-size scaling of ^(Tr, 7r/2) for U/t = 9.0 and 
t'/t = 0.7 at half filling. The result suggests 2 x 4 AF long- 
range order. Open circles for 6x6 and 10 x 10 lattices are 
the peak values while filled symbols are the intensity of peak 
summed over k points at (7r,7r/2) and its nearest neighbor k 
points. 



model in the leading order by 

H — Ji'Ej^ij-^Si ■ Sj + J2^{{ij))Si ■ Sj (8) 

where () and (()) denote the nearest- neighbor sites and 
next-nearest-neighbor sites, respectively. Here, Ji = 
At'^/U and J2 = 4i'^/[/ are both antiferromagnetic 
interactions so that magnetic frustration arises. For 
J2/J1 < 0.4 the Neel order with peak structure in S{q) 
at g = (tt, tt) was proposed, which corresponds to the 
AFI phase in the Hubbard model. On the other hand, 
for J2/J1 > 0.6, the stripe order with q — (0,7r) or 
(tt, 0) peak in S{q) was proposed, which corresponds to 
the stripe shape jsj in the next section. For the inter- 
mediate region of 0.4 < J2/J1 < 0.6, no definite con- 
clusion has, however, been drawn on the nature of the 
ground state. The possibihty of the columnar-dimerized 
state [23, I23, I23, the plaquette singlet state |3Q| and 
the resonating-valence-bond state were discussed. In 
the present study, we find new type of AFI phase with 
longer period between AFI and stripe phases. This new 
AFI phase could have some connection to the region 
0.4 < J2/J1 < 0.6. In terms of the effects of frustra- 
tions, the stabilization of the longer-period AF structure 
is a natural consequence in the classical picture, where 
well known ANNNI model exhibits a devil's staircase 
structure [SJ. Complicated longer-period structure may 
melt when quantum fluctuations are switched on, while 
the present result indicates the survival of 2 x 4 structure 
for large U. 

To investigate the possibility of the dimerized state and 
the plaquette state, we consider the dimer correlation 



where 



Oa 



1 ^ 



i-}-a 



(9) 



(10) 



i=l 



yy 



and a shows the unit vector in the a direction. D 
clearly indicates that the dimer order in the y direction 
is absent as we see in Fig. [7\ On the other hand, in this 
definition of D, D^x should also show the long-range or- 
der, if the 2x4 AF order is stabilized. This is indeed 
seen in our data. In the limit of strong coupling, this 
region is mapped to the Heisenberg model with nearest- 
neighbor exchange Ji and the next-nearest-neighbor ex- 
change J2, with J2/J1 ~ 0.5. In this region, columnar 
dimer state has been proposed as the candidate of the 
ground state [23, 123, 123 • Our result also shows the dimer 
long-range order supporting the existence of the colum- 
nar dimer phase. However, in this phase at finite U, 
the antiferromagnetic long-range order with longer pe- 
riod with 2x4 structure shown in Fig. 131(B) is also seen. 
Since the dimer-order parameter defined by D^x is au- 
tomatically nonzero in this extended antiferromagnetic 
phase, the primary order parameter should be identified 
as the longer-period antiferromagnetic order. We nat- 
urally expect a continuous connection of the Hubbard 
model to the Heisenberg model, while as far as we know, 
serious examination of such longer period AF order is not 
found in the literature. It is an intriguing issue to exam- 
ine this new possibility of longer-period AF order in the 
Heisenberg limit. 



VII. STRIPE PHASE 

For larger t'/t and U/t, another phase appears, whose 
basic configuration is stripe shape in Fig.|3|[C). In Fig. |S1 
we show strong AF bonds for t and t' directions. For AFI 
phase, spins on the bonds along t-direction become anti- 
ferromagnetic each other and gain energy while spins on 
the bonds in t'-direction become parallel and lose energy. 
On the other hand, for the stripe shape, spins on the 
bonds for i'-direction become all anti-parallel with com- 
promised anti-parallel spins on the bonds for i-direction. 
In this situation, as t'/t becomes larger, stripe phase be- 
comes energetically favored. As shape (B) is between 
(A) and (C), it becomes energetically favored at medium 
value of t'/t. 

The minimum block of stripe shape is 2 x 1. Therefore, 
there is no irregularity of spin correlation for different 
lattice sizes. The spin correlation has a sharp peak at 
(0,7r) as we see in Fig. EJC). Finite-size scaling in Fig|3 
indeed shows the presence of the long-range order in the 
thermodynamic limit. The overlap between basis state 
and its shifted one also indicates that the configuration 
in Fig. |3IC) is realized. 
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FIG. 9: Finite-size scaling of 5(0, tt) for f//f = 9.0 and f'/f = 
1.0 at half filling. The result suggests that the stripe has 
long-range order at this parameter value. 



FIG. 7: (color online) System-size dependence of the dimer- 
correlation functions (Dx:c and Dyy) for t = 1.0, t'= 0.7 and 
U = 9.0 at half filling on iV = 4 x 4, 6 x 6, 8 x 8, 10 x 10 and 
12 X 12 lattices. The blue open circles, the red open diamonds 
and green open squares show the D:cx and Dyy of AF2, Stripe 
and AF phases, respectively. The long-range order suggested 
in Dyy implies the 2x4 AF order. 
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FIG. 8: (color online) Strong AF bonds of AFI phase (A), 
new AFI phase (B) and stripe phase (C). 



VIII. NON-MAGNETIC INSULATOR PHASE 



Three phases of previous sections are semi-classical. 
Their main configurations can be intuitively understood 
from the classical picture of the spin alignment. Next 
we consider an unconventional phase of frustrated Hub- 
bard models, which is illustrated in Fig|2] as the non- 
magnetic insulator (NMI). This phase does not have any 
distinct peak structure in 5(k) as we see in Fig.|SjD). Be- 
cause this NMI phase appears in a window sandwiched 
by the metal and insulating antiferromagnetic phases, it 
is clearer than the previous studies 5, 6] that the phase is 
stabilized by charge fluctuations enhanced near the Mott 
transition. 



A. spin excitation gap 

Now in the NMI phase, system size dependences of 
the spin excitation gap A_E between the ground state 
and the lowest triplet are shown in Fig. ^1 Figure [Till 
indicates that the triplet excitations become gapless in 
the thermodynamic limit. The gap appears to be scaled 
asymptotically with the inverse system size N~^, namely 
AE ^ C/N. At least the extrapolated gap (< O.Oli) is 
much smaller than the typical gap inferred in the spin 
gapped phase of the corresponding Heisenberg limit (> 
O.IJ) |2J|. The gapless feature shares some similarity to 
the behavior in the AFI phase. However, the detailed 
comparison clarifies a crucial difference as we will show 
later. The fitting in the NMI phase to the form (6) gives 
unphysical values such as c > 1.5. We note that the 
uniform magnetic susceptibility is given by 2/3C, which 
implies a nonzero and finite uniform susceptibility. 

Except for ID systems, the present result is the first 
unbiased numerical evidence for the existence of gap- 
less excitations without apparent long-ranged order in 
the Mott insulator. Although a tiny order cannot be 
excluded if it is beyond our numerical accuracy, in the 
present NMI phase, the absence of various symmetry 
breakings including the AF order has already been sug- 
gested by the size scalings of the models [A] jig and 
[B] 6,]. It was shown 25] for the model [A] that as well 
as dimer and plaquette singlet orders, s- and d-density 
waves are also numerically shown to be unlikely for four 
types of correlations probed by 

CM) = |(J«(q)4(q))| 

^(q) = ^E'^L'^k+q,./.(k) (11) 
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FIG. 10: (color online) Size scalings of the 5 = 1 excitation 
gaps for several choices of parameters. The triangles show the 
case of the model [B] while others are for the model [A] . The 
dashed curves are fittings to Eq.(6). The circles, squares and 
triangles are results obtained inside the NMI phase. 



with 



/i(k) = cos(fcj;) + cos(/cj,), 

/2(k) = cos(fc^) - cos(fcj^), 

/3(k) = 2cOs(fca;)cOs(fcy), 

/4(k) = 2sin(fc2;)sin(fcy). 



(12) 



The calculated results all indicate the absence of long- 
ranged orders with the correlations staying rather short 
ranged. However, it does not strictly exclude an ex- 
tremely small and nonzero order parameter, although it 
may melt with further decreasing [/. 



B. Lowrest 5 = 1 and 5 = excitations 

The lowest energy S = \ excitations at each momen- 
tum sectors -Ei(k) show further dramatic difference be- 
tween the AFI and NMI phases. The lowest energy states 
with specified momenta, i?(k) are calculated from the 
spin-momentum resolved PIRG. When £'(k) becomes the 
maximum at ^max and the minimum at k„ii„, we intro- 
duce the width W = E{\s.^ax) - ^(k„i„). 

We first make a general remark on W expected from 
the known phases. Note that in the AFI phases as well 
as in the stripe phase, the lowest S' = 1 excitations i?i(k) 



give nothing but the spin wave and the momentum de- 
pendence is given by the spin wave dispersion. There- 
fore, W represents the dispersion width. On the other 
hand, in the Fermi liquid, it is given by the lowest edge 
of the continuum of the Stoner excitations (particle-hole 
excitations) and in a certain range of the momenta, they 
are gapless. For example, for the case of noninteracting 
Fermions at half filling on the square lattice, W becomes 
zero. In numerical calculations of the Fermi liquid, how- 
ever, it has, in general, a finite width W > ^ due to the 
finite-size effect and the width is expected to vanish as 
W oc XjN^I'^ ^ namely being scaled by the inverse linear 
dimension of the system size. 

In the AFI phase, the dispersion is essentially described 
by the spin-wave spectrum similar to 



l^E(k) ^AJJl 



ll 



with 



7fe = 2 (cos(fci) + cos(A:j^)) 



(13) 



(14) 



for the spin wave theory of the Heisenberg model, but 
modified because of finite U . The calculated dispersion 
width at [/ = 4,t = 1, and t' = 0, shows a small system 
size dependence and is around 1.5, which can be com- 
pared with the dispersion width (^ 1.6 for J = 0.4) of 
the ordinary spin wave obtained from the map ping of the 
Hubbard model to the Heisenberg model [2m |. 

In marked contrast, the width W for 5 = 1 excita- 
tion in the NMI phase has strong and monotonic system 
size dependence as in Fig. 1111 For systems larger than 
8x8 lattice, the dispersion becomes vanishingly small. 
The vanishing W is consistent with what expected in 
the Fermi liquid, although the spin liquid phase is cer- 
tainly insulating. It implies that only the "spinon" ex- 
citations form an excitation continuum in the presence 
of the charge gap. Although it is not definitely clear, 
the size dependence seems to show very quick collapse 
of the dispersion with increasing system size and may 
not be fitted by a power of the inverse system size as 
in the single-particle Stoner excitations in metals. Such 
quick collapses of W are observed solely in the NMI 
phase irrespective of the models (namely commonly seen 
in the models [A] and [B]). The collapse implies that 
the triplet excitations cannot propagate as a collective 
mode. We also note that the gap of S' = 1 excitations 
from the ground state, namely AE is scaled by 1/N as 
described above. Therefore, the momentum degeneracy 
within S = 1 sectors seems to be much higher than the 
spin degeneracy in the thermodynamic limit. 

The presence of such degenerate excitations well ac- 
counts for the quantum melting of simple translational 
symmetry breakings including the AF order, because 
a long-ranged order in the two-dimensional systems 
is destroyed when the excitation becomes flatter than 
AE{k) = k"^ . This is because of the infrared divergence 
of fluctuations in the form / k'^~^dk-^^T^ with d = 2. 
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FIG. 11: (color online) Size scalings of W for S = and S = 1 
excitations in the NMI phase for the model [A]. 



The total singlet state {S = 0) at any total momentum 
k also shows degenerate structure in the ground state for 
larger system sizes as in Fig. ^J We similarly introduce 
the width W as W = ECkmax) — E{V.min) within the 
singlet S" = sector. The width W vanishes in the NMI 
phase in the model [A] as well as in [B] . This again implies 
the excitation continuum, which has larger degeneracy 
than the spin excitations. 



C. Spin renormalization factor 

The spin renormalization factor Zg{q) is defined as 

Zs{q) = \{S^l,q\Z\S = Q,q = 0)\ (15) 

where 

^ = ^ II (4+g,TCfc,T - 4+g,i^^a) • (16) 

It is rewritten in the site representation as 



<"-71 



J 



{S=l,q\{r 



-Tj 



Hj, 



15 = 0,9 = 0), 



-iqRj 



(17) 
The PIRG wave functions for 5 = and 5 = 1 are 
given by spin-projection operator. As the derivation of 
its matrix element is somewhat lengthy and needs spin 
algebra, it is summarized in Appendix. 

In analogy with the renormalization factor of the quasi- 
particle weight in the Fermi liquid, Zs measures whether 
the spin excitation is spatially extended and can propa- 
gate coherently or not. If Zs has nonzero values, the spin 
excitation can propagate coherently. Figure 1121 indeed 
shows that Z^ remains nonzero for these ordered phases. 



In fact, in AFI, AFI2 and stripe phases, the magnons are 
well-defined elementary excitations in momentum space 
and spatially propagates coherently, which are reflected 
in nonzero values of the extrapolated Zs to the thermo- 
dynamic limit. In contrast, the renormalization factor 
appears to scale to zero for NMI phase, which implies 
that the spin excitations dressed by other spins are spa- 
tially localized. 
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FIG. 12: (color online) Size dependence of spin renormaliza- 
tion factor and S{q) at the peak value are shown for AFI, new 
AFI(AFI2), Stripe and NMI phases. 



IX. DISCUSSIONS AND SUMMARY 

The present excitation spectra show the following 
double-hierarchy structure: Although the singlet ground 
state is unique, there exist enormous number of low- 
energy excitations within the 5 = sector leading to 
nearly dispersionless momentum dependence of singlet 
excitations. Although the energy of lowest-energy state 
at each total momentum has a substantial momentum de- 
pendence in finite-size systems, thus forms a dispersion- 
like structure arising from the finite-size effects, the mo- 
mentum dependence quickly collapses with increasing 
system size. Another degeneracy near the ground state 
appears in the excitations with 5 > 0. The lowest-energy 
states with 5 = 1, 5 = 2, • • • denoted hy Ei,E2,- ■ ■ have 
a tower-like structure, Ei < E2 < ■■■ in finite-size sys- 
tems. However, these spin excitation gaps again collapse 



with increasing system size supporting gapless spin ex- 
citations in the thermodynamic hmit. With increasing 
system sizes, the singlet excitation energies become van- 
ishing with much faster rate than that of vanishing exci- 
tation energies for the excitations to nonzero spins such 
as triplet. In other words, the collapse of the excita- 
tion energies for the spin dependence seems to be much 
slower than the collapse in the momentum dependence. 
This generates a hierarchy structure of the excitation en- 
ergies. 

Here we discuss a possible interpretation of the prop- 
erties. The vanishing gap may allow the following inter- 
pretation: The ground state is given by a superposition 
of dynamical singlet bonds, which cover the whole lat- 
tice. Those singlet bonds with vanishingly small singlet 
binding energy have a nonzero weight, because of the 
distribution of the singlets over long distance, where, for 
example, weights decays with a power law with increas- 
ing bond distance as in a variational long-ranged RVB 
wavef unction j33 |. 

The collapse of W and vanishing renormalization fac- 
tor Zs show that the ground state is nearly degenerate 
with other low-lying states obtained by spatial transla- 
tions and they have vanishing off-diagonal Hamiltonian- 
matrix elements each other. Although the translational 
symmetry is retained in the original Hamiltonian, the 
present orthogonality may imply a large degeneracy near 
the ground state. The absence of the spin renormal- 
ization factor Zs = implies first that the Goldstone 
mode like magnons in the symmetry broken phase does 
not exist. Furthermore, it also excludes the existence 
of coherent single spin excitations as in the idea of the 
particle- hole- type excitations at the hypothetical spinon 
Fermi surface |35ll3q|. The result rather suggests that a 
spin excitation obtained from the unbinding of the weak 
singlet may be spatially localized because of the dressing 
by the surrounding sea of singlets. 

In classical frustrated systems, the spin glass can be 
stabilized by an infinitesimally small quenched random- 
ness, for example, in the Ising model on a triangular lat- 
tice |33, where the macroscopic degeneracy remains in 
the regular system. In the present case, it may also be 
true that tiny randomness may further stabilize the freez- 
ing of the localization of spins and leads to the spin glass. 

The nature of the gapless and dispersionless excitations 
is not completely clarified for the moment. Although the 
coherence of the spin excitations must be more carefully 
examined, the present result supports that an unbound 
spin triplet does not propagate coherently due to strong 
scattering by other weakly bound RVB singlets. This 
result also means that the quantum melting of the AF 
order occurs through the divergence of the magnon (or 
Goldstone mode) mass. 

This NMI phase appears to be stabilized by the Umk- 
lapp scattering [3»|, where the Mott gap is generated 
without any symmetry breaking such as antiferromag- 
netic order. The degenerate excitations within the sin- 
glet, which are similar to the present results, but with a 



spin gap were proposed in the Kagome and pyrochlore 
lattices based on small cluster studies j39|]. The possible 
symmetry breaking from degenerate sing lets was also ex- 
amined on the pyrochlore lattice [i^, |41| , while the spins 
were again argued to be gapful in contrast to the present 
results. In our results, the gapless spin excitations be- 
come clear only at larger system sizes than the tractable 
sizes of the diagonalization studies. 

We briefly discuss experimental implications of the 
present new quantum phase. Recent results by Shimizu 
et al. [i3 on k-(ET)2Cu2(CN)3 appear to show an ex- 
perimental realization of the quantum phase we have dis- 
cussed in the present work. In fact this compound can 
be modeled by a single band Hubbard model on nearly 
right triangular lattice near the Mott transition. The 
NMR relaxation rate shows the nonmagnetic and gapless 
nature retained even at low temperatures (^ 0.2K) and 
suggests the present quantum phase category. Another 
organic compound also shows a similar behavior |ll| 

Systems with Kagome- (or triangular-) like 
structures, '^He on graphite [13 and volborthite 
Cu3V207(OH)2-2H20 |T3 show nonmagnetic and gap- 
less behaviors. On the other hand, glass-like transitions 
are seen in 3D systems, typically in pyrochlore com- 
pounds as i?2M0207 with R =Er, Ha^Y, Dy, and Tb ^ 
and in fee structure, Sr2CaRe06 [3- It is remarkable 
that the glass behavior appears to occur even when ran- 
domness appears to be nominally absent in good quality 
single crystals of stoichiometric compounds. Although 
the lattice structure, dimensionality and local moments 
have a diversity, many frustrated magnets show gapless 
and incoherent (glassy) behavior. The present result 
on gapless and degenerate structure emerging without 
quenched randomness offers a consistent concept with 
this universal trends. It would be an interesting open 
issue whether the present system leads to such a glass 
phase at T = in 2D by introducing randomness. 

In summary, we have studied the 2D Hubbard model 
on two types of square lattices with geometrical frustra- 
tion effects arising from the next-nearest neighbor trans- 
fer t' . In the parameter space of the interaction U and 
t', paramagnetic metal, two antiferromagnetic insula- 
tor phases, a stripe-ordered insulator phase, and a new 
degenerate quantum spin-liquid phase are found. The 
quantum spin liquid (in other words, nonmagnetic insu- 
lator) phase has gapless spin excitations from the degen- 
erate ground states, and furthermore the dispersionless 
modes are found in all the spin sectors. The calculated 
spin renormalization factor suggests that the gapless spin 
excitation is spatially localized and does not propagate 
coherently in the thermodynamic limit. Recent exper- 
imental findings, including quantum spin liquids in 2D 
and spin glasses in 3D suggest a relevance of this phase 
in disorder-free and frustrated systems. 
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Appendix 

In this appendix, we discuss how to evaluate the ma- 
trix elements between different spin states. We con- 



sider the state with definite spin S and its z-component, 
which is given from general Slater determinant |(/)) by 
spin projection Ps,m- The detailed form of Ps,m is given 
by J2k 9kLm ji in Ref.|23, HI] where gx's are param- 
eters. We denote such spin projected wave function as 
\(J)s,m) = Ps,m\4')- Now we consider an operator O^ 
which changes spin quantum number by A and its z- 
components by ^. Its expectation value is evaluated by 
Wigner-Eckart theorem. Therefore we introduce reduced 
matrix element (a, S \\0^\\ a' , S") generally defined by 



J 



{a,S,M\0^^\a',S\M') = {-) 



S-M 



S X S' 
-M n M' 



, {a,S\\0^\\a',S') 



(18) 



where a and a' are additional quantum numbers. By considering the commutation relation between spin-projector 
and the operator Ofj, a following formula can be derived |42| as 



/,n lln^ll w!, \ _ (2Si + l)(2So + l) \^ „* „ ( \So-Ko f So A Si 

ifs, \\0 II (pso) 8^F2 2^ 9ko9k^ {-) I _^^ ^ ^^ 



KoKiKi,, 



dnD'=.\^(^\o':R{n)\^) (i9) 



where, in reduced matrix element, the z-component is suppressed, Dfj ^ is Wigner's _D-function and R (fl) is a 
rotation operator in spin space p^ . 

Now we consider the spin renormalization factor. Its operator O}^ is defined by 



oi = E(i^'i'^'|i'^)4c. 



(20) 



where Cj^^r = {—)^ '^Cj-^. Its three components are Oq = -y= I — clc| -I- cjcj^j, O} = clc^ and O^i — — c|c|. As the 
wave function is determined in the half-filled space, the z-component of spin is zero, which means gx = Sk,q- Moreover 



we set 5o = and 5*1 = 1 and normalization factor of projected wave functions like y - — j^' ' - — ^^-^ ^ — 
taken into account. Therefore we can obtain a following formula as 

((^1,0 \Ol\ 0o,o) = E h^^PdMU iP) iv \0y^''-\ 4>)- 
In the nuclear structure physics, the formula Ea. fTi^ is ofte n used. 



^ is also 



(21) 
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